// random number generators using Mersenne twister and polar Box-Muller transform
// algorithms taken from en.wikipedia.org/wiki/Mersenne_twister and en.wikipedia.org/wiki/Box-Muller_transform
#include <ClanLib/core.h>
#include "random.h"
#include <time.h>
#include <cmath>

Random::Random()
{
	index = 0;
	count = 0;
	// use current time as seed
	this->seed((unsigned int)time( NULL ));
}

Random::Random(unsigned int seed)
{
	index = 0;
	count = 0;
	this->seed(seed);
}

Random::~Random()
{
}

// seed the random number generator
void Random::seed()
{
	// use current time as seed
	this->seed((unsigned int)time( NULL ));
}

void Random::seed(unsigned int seed)
{
	//std::cout << "Seed random numbers" << std::endl;
	MT[0] = seed;
	for(unsigned short i = 1; i < 624; i++)
	{
		MT[i] = (1812433253ULL /* 0x6c078965 */ * (MT[i-1] ^ ((MT[i-1]) << 30)) + i) % 4294967296ULL; // 0x100000000
	}
}

// generate array of 624 untempered numbers
void Random::refreshNumbers()
{
	//std::cout << "Refresh random numbers" << std::endl;
	unsigned int num;
	for(unsigned short i = 0; i < 624; i++)
	{
		num = ((MT[i])&(1 << 31)) + ((MT[(i+1) % 624])&~(1 << 31));
		MT[i] = MT[(i + 397) % 624] ^ (num >> 1);
		if((num % 2) == 1)  // num is odd
		{
			MT[i] = MT[i] ^ (2567483615); // 0x9908b0df
		}
	}
}

// generate uniform random unsigned integers between 0 and RANDOM_MAXRAND
unsigned int Random::rand()
{
	if(index == 0){
		refreshNumbers();
	}

	int num = MT[index];
	num = num ^ (num << 11);
	num = num ^ ((num >> 7) & (2636928640)); // 0x9d2c5680
	num = num ^ ((num >> 15) & (4022730752)); // 0xefc60000
	num = num ^ (num << 18);

	index = (index + 1) % 624;
	return num;
}

// generate uniform random unsigned integers between 0 and max
unsigned int Random::rand(unsigned int max)
{
	unsigned int result = (unsigned int)(((double)rand() / RANDOM_MAXRAND) * (max + 1));
	if(result > max) result = max;
	return result;
}
// generate uniform random unsigned integers between min and max
unsigned int Random::rand(unsigned int min, unsigned int max)
{
	unsigned int result = (unsigned int)(rand() * (max - min + 1) / RANDOM_MAXRAND) + min;
	if(result > max) result = max;
	return result;
}

// generate uniform random numbers between 0.0 and 1.0
double Random::uniform()
{
	return (double)rand() / RANDOM_MAXRAND;
}

// generate uniform random numbers between 0.0 and max
double Random::uniform(double max)
{
	return max * uniform();
}

// generate uniform random numbers between min and max
double Random::uniform(double min, double max)
{
	return min + (max - min) * uniform();
}

// generate gaussian random number with mean 0 and standard deviation 1
double Random::gaussRand()
{
	if(gaussStored)
	{
		gaussStored = false;
		return tempGauss;
	}

	double u, v, s, multiplier;
	do {
		u = uniform(-1.0, 1.0);
		v = uniform(-1.0, 1.0);
		s = u*u + v*v;
	} while(s >= 1 || s ==0);
	multiplier = sqrt((-2 * log(s)) / s);
	tempGauss = (v * multiplier);
	gaussStored = true;
	return (u * multiplier);
}

// generate gaussian random number with mean 0 and given standard deviation
double Random::gaussRand(double stDev)
{
	return stDev * gaussRand();
}

// generate gaussian random number with given mean and standard deviation
double Random::gaussRand(double mean, double stDev)
{
	return mean + (stDev * gaussRand());
}
